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ABSTRACT 

The Palomar High-precision Astrometric Search for Exoplanet Systems (PHASES) has monitored 
37 sub-arcsecond binary systems to determine whether their Keplerian orbits are perturbed by faint 
astrometric companions to either star. Software has been developed to evaluate the regions in a 
companion mass-period phase space in which the PHASES observations can exclude the possibility of 
face-on orbit perturbations. We present results for 8 systems for which astrometric companions with 
masses as small as those of giant planets can be excluded. 

Subject headings: binaries: close - binaries: visual - planetary systems: formation - astrometry - meth- 
ods: data analysis 
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Q\ • 1. INTRODUCTION 



A technique has b een developed to obtain high p recision (10-20 /zas) astrometry of close stellar pairs (separation less 
than one arcsecond: lLane and Muterspaug h 2004) using long-baseline infrared interferometry at the Palomar Testbed 
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Interferometer (PTI: 



Colavita et alT l999). These observations provide precision visual orbits of the binaries and allow 



■ detection of tertiary components orbiting either the primary or secondary due to the reflex motion of the subsystem 
' center-of- light. 

, These measurements were made at PTI as part of the Palomar High-precision Astrometric Search for Exoplanet 
qq " Systems (PHASES) program. PTI is located on Palomar Mountain near San Diego, CA. It was developed by the 
Jet Propulsion Laboratory, California Institute of Technology for NASA, as a testbed for interferometric techniques 
\Q , applicable to the Keck Interferometer and other missions such as the Space Interferometry Mission (SIM) . It operates 

■ in the J (1.2/im), H (1.6/im), and K (2.2/im) bands, and combines starlight from two out of three available 40-cm 
apertures. The apertures form a triangle with one 110 and two 87 meter baselines. PHASES observations use two of 
the three available baselines at PTI: the NS (110 meter) and SW (87 meter) baselines. 

PTI has previously contributed to exo planet search efforts through high spatial resolution visibility "imaging" of 
the 51 Pegasi svstem (|Boden et al.lll998T) . which placed limits on the luminosity of its companion and support the 
+3 . conclusion that it is substellar, and likely planetary, in nature. Here the initial detection limits of a different effort at 
52 ■ PTI are presented, this time based on astrometric measurements of binary systems. 

2. THE PHASES EXOPLANET SEARCH 

PHASES measures the separation vectors of bright binaries that are not resolved by PTFs telescopes but are over- 
rS | resolved by the interferometer itself. The high spatial resolution (~ 4 mas) of the interferometer, coupled with an 
hour of observations, enables extremely precise measurement of the binary separation. Phase-referencing is used to 
maintain coherence so that the full resolving power of the interferometer is applied to the astrometric measurement. 
The current limiting magnitude in phase-referencing mode is K = 4.5 for equal magnitude binaries. Work is being 
done to improve the sensitivity by introducing a variation in which the phase-referencing fringe tracker operates at 
half the current speed, to a 50 Hz duty cycle. The PHASES program has successfully obtained at least one night of 
data on 37 binaries. As of April 2006, observations have been attempted on 209 nights, resulting in 688 astrometric 
measurements . 

If only the primary and secondary stars are present in the system, one expects these separation vectors to evolve 
according to a Keplerian model. On the other hand, if additional (faint) components are present, their presence will 
cause perturbations to this orbit. Distant companions simultaneously orbiting both visible components will cause only 
very small perturbations to the observed separation vector (the differential gravitational pull is small); the PHASES 
obser vations are not sensitive to these "circumbinary" planets (also called "P-type" or planetary-type planets, iDvorakl 
However, companions orbiting just one star of the binary can cause noti ceable perturba tions to the separation 
vectors. This configuration is similar to that outlined by the Nemesis theory l|Harrisonlll977j) . which postulated the 
Sun has a stellar companion orbiting at a distance far from the planets. PHASES is a search for these "Nemesis-type" 
planetary systems (also called "S-type" or satellite- type planets). 

The primary goal of the PHASES program is to find and characterize giant planets in close binary systems. In this 
case, "close" binary is defined as those with semimajor axis a < 50AU in which one might expect binary dynamics to 
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play a major role in system formation and evolution (see, for example, Pf ahl and Muterspaug h 2006). The existence 
of such systems poses a strong challenge for models of gia nt planet formation. While it is possible each of the two 
processes currently favored — core accretion l)Lissauerlll993|) and gravitational instability l)Bossll2000|) — contribute to 
giant planet formation around single stars and wide binaries, simulations show both sch emes have ob stacles when a 
second star orbits so closely that it interacts with the planet-forming circumstellar d isk ( Nelson 2000; ), truncating i t 
in size and heating it . In five exoplane t hosting bina ries — HD 188753 llKonackill2005D. 7 Cephei l)Hatzes et alJ l20()3). 
GJ 86 l|Queloz et al.ll2000j) . HD 41004 l|Zucker et aUl200|) . and HD 196885 (|Chauvin et al. 112006]) — the secondary star 
would have truncated the disks to less than 7 AU (for systems in which the binary orbits are not fully constrained, 
moderate values for the eccentricities are assumed) . It is possible that some of these planets formed around single stars 
or in wide binaries and reached their current configurations via dynam ical interactions in the short-lived star clusters in 
which they formed (Pfahl 2005; Portcgics Zwart and McMillan 2005), though this post-formation mechanism appears 
to be too infrequent to explain the number found. 

These systems have been identified with the radial velocity (RV) method. PHASES employs astrometry, from which 
the companion mass can be identified without the ambiguity of the orbital inclination. Furthermore, the relative 
orientations of the binary and planetary orbits can be determined in order that system dynamics can be studied. This 
effort and others specifically targeting close binaries will better determine the frequency of planets in binaries; if large, 
this will be strong motivation for a revolution in the theory of giant planet formation. 

3. SYSTEMATIC ERRORS 

When fit to a low-order polynomial or two-body (single) Keplerian model, the current PHASES measurements show 
night-to-night scatter in excess of that determined by evaluating the scatter within a given night. The timescales of 
this excess scatter can be evaluated by dividing observations from within a night into subsets by time and analyzing 
these subsets individually. It is found that within a given night the subsets agree at the level determined by the 
formal uncertainties obtained from the standard PHASES data reduction pipeline. Because the subsets represent 
observations separated by of order an hour, it is concluded that the PHASES uncertainties are consistent within a 
night; the variations occur on night-to-night timescales. 

This excess scatter may be due to either the presence of additional companions or due to systematic errors that 
occur on timescales of days or longer. In order to determine limits to tertiary companions, it is assumed that any 
systematic errors result in true night-to-night measurement uncertainties that are related to the formal uncertainties 
by a multiplicative factor. Depending on the nature of potential systematics, this factor may be different for each 
individual system (potentially depending on contrast ratios, separations, or star color differences); it is assumed that 
this factor is a constant for any given system over the course of the measurements. An alternative approach would 
be to postulate a noise floor; this is not evaluated in the current investigation. It is assumed the uncertainties are 
distributed in a Gaussian manner; see section l5~4l 

Potential sources of systematic instrumental errors have been identified and efforts are being made to correct them 
for the next stage of PHASES observations (during the 2006-2008 seasons). The first source is due to color dispersion 
within the interferometer, which manifests itself in two manners. The path compensation for the geometric delay at 
PTI has been done with movable mirrors in air, which has a wavelength-dependent index of refraction. The fringe 
packets of astrophysical sources are dispersed by an amount that depends on the difference in air paths between arms 
of the interferometer; this changes the shape and overall location of the fringe packets. If the two stars are of differing 
colors, each will be dispersed by a slightly different amount, and their apparent separation will be biased; see the 
schematic in figure ^ Additionally, if there is dispersion in a direction perpendicular (lateral) to the light beam (due 
to optics acting as prisms; atmospheric differential chromatic refraction also introduces lateral dispersion, though this 
amount is small) or diffraction, the color of the star's light falling on the image plane detector pixel can vary with sky 
position. In the presence of longitudinal dispersion, the dependence of these color shifts on location within the pixel can 
also introduce astrometric systematics; see figure [5] A longitudinal dispersion compensator has been commissioned for 
PTI which should eliminate both systematics. A second source of potential error is drift in the "astrometric baseline" 
due to drifts in optical alignment which may result in variable pupil sampling at the interferometer apertures. This 
is remedied by introduction of an automatic alignment system. These engineering improvements will initiate a second 
stage of the PHASES program, in which it is anticipated the observational precisions will be improved by a factor of 
~ 3. 

4. DETECTION LIMITS 

It is anticipated that at least a portion of the excess noise is the result of systematic errors, though in some cases 
the presence of an astrometric companion also contributes. An approach to differentiate between these contributions 
is developed by recognizing that an astrometric perturbation will appear as a Keplerian wobble, whereas instrumental 
terms will be random; the Keplerian signature becomes a constraint on the nature of the excess scatter. An algorithm 
for determining the confidence levels of such a signal in the presence of instrumental scatter is described below. 

4.1. Exclusion Criteria 

The hypothesis is that model A correctly represents the observed system; one desires to determine the \ 2 goodness- 
of-fit threshold at which a different model, B, is less consistent with the data than A (for a discussion of x 2 as a 
likelihood estimator for data with Gaussian uncertainties, see for example iPress eTaDll99^ . In the present case, 
model A is the simpler single-Keplerian orbit (with seven independent parameters allowed to vary freely, P\, ex, i\, u>i, 
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Fig. 1. — Schematic of the shift in fringe positions due to dispersion (the effect has been exaggerated for clarity). The vacuum (no 
dispersion) interferograms are plotted with solid lines; those dispersed by air with dotted lines. 

(top) Dispersion shifts the point of zero optical path difference for a star, due to different amounts of air path in each arm of the 
interferometer (the effective optical path difference measured as if in vacuum). 

(middle) The dispersion shifts for stars of equal colors are equal and cancels; the measured separation is the same. 

(bottom) Stars of unequal colors are shifted by slightly different amounts by dispersion, and the resulting measured separation is different. 
Not shown are the shape distortions to interferograms. In a differential measurement, these cancel to first order and are insignificant at 
the precisions of PHASES measurements. 
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Fig. 2. — Schematic of the shift in fringe positions due to coupling of longitudinal and lateral dispersion. Due to the finite distance 
between stars in the image plane, the detector will sample each star differently. If the stars' spectra are dispersed in the image plane, 
a given detector pixel will sample different colors for each star. This is normally not a problem, as the astrometric measurement is not 
derived from centroiding the stars on the detector, but rather from the locations of the movable mirrors at which fringes appear on the 
detector from either one star or the other (each star image is the overlap of images from two telescope). If there is additionally longitudinal 
dispersion, the delay location at which fringes form will be color dependent. 



T a i, 0,%, and a\ ), whereas model B is the double-Keplerian model with additional parameters. In the case of a second 
"perturbation" orbit being limited to face-on circular configurations, three new free parameters will be introduced: 
period P2, center-of-light semimajor axis ai 1 and one member of the set of the angle of periastron passage 0J2, epoch 
T 0i 2, and longitude of the ascending node O2 , which are degenerate; constant parameters are the other two of the 
previous set, the inclination 12, which can either take value zero or 1 80 d egrees (though nothing in between), and 
the eccentricity e2 = 0. Extension to inclined orbits is discussed in A phase space can be constructed over 

these additional free parameters {P2, 02, and O2 ) to evaluate which combinations result in improved or worse fits. 
In practice, this phase space is collapsed in the il 2 dimension to the best-fit value, and 02 is converted to tertiary 
companion mass. 

Qa and Qb are defined as the number of free parameters for models A and B, respectively (here, Qa < Qb), and 
there are D independent measurements being considered (for astrometry, D is twice the number of measures, as each 
are two dimensional). Thus, the numbers of degrees of freedom are Za = D — Qa and Z& = D — Qjg. 

In x 2 fitting, the best-fit model is expected to have a goodness-of-fit \ 2 equal to the number of degrees of freedom. 



4 



In the hypothesis that model A accurately describes the data, but that the data have excessive noise not included 
in their formal uncertainties, the noise excess factor is given by Xa/Za- Reweighting the formal uncertainties by the 
square root of this factor will result in a recomputed value of Xa mod = Za- The hypothesis thus assumes that this 
new set of uncertainties correctly represents the scatter in the observations as a random process. 

When introducing additional free parameters to a model, it is expected that the fit x 2 will decrease. This does not 
mean that the model fits better, rather this must be placed in terms of the reduced x 2 = X 2 /Z- If the two models 
represent the data equally well, one expects 

Xa/Za = xI/Zb- 

If this relationship does not hold and B is has the larger reduced x?) B is excluded at the Na level for the value of 
N that satisfies 

X 2 a/Za = XbI ( z b + Sq b ,n) 

where <Jq,jv is the value by which a properly normalized \ 2 metric must be increased to find the Na confidence region, 
given a model with Q free parameters. In this case, it is the interval computed for model B, which by assumption less 
accurately describes the data. 
Rearrangement of terms indicates that model B is excluded at or beyond the Na level when 

Xl > Zb+ z Sq - N X 2 a- (1) 
"A 

4.2. Astrometric Detection Criteria 

In this hypothesis, we accept B as the correct model, and look for the \\ level for which model A is excluded. Thus, 
we simply replace the roles of the models in equation^ 

2 ^ Z A + Sq a . N 2 

Xa 7 Xb- 

^B 

Because it is model B for which a phase space grid is constructed (over the extra free parameters in this more 
complex model), it is useful to invert this expression such that it is expressed as the contours of x% f°r which model 
A is rejected: 

xl < tt4j Xa- (2) 

Tabled lists values of 8q,n used in the following sections. The values of (5q,at can be found by iteratively solving 
the equation 

P (Q^) =P (l ^ 



where 

P(a,x) 



2 ' 2 J V 2 2 
J x e-H^dt 



is the standard incomplete gamma function, equal to the probability of a \ 2 distribution. In practice, it is often 
numerically better to equate Q(a, x) = 1 — P(a, x). Thus, the la contour has probability 0.683, and 3<r has probability 
0.9973. 

One must be careful how the exclusion and detection confidence levels are interpreted. This is particularly true 
in the case of detections; one must not forget the effect of the data sampling function or the possibility that a third 
model better explains the data. To claim a true astrophysical "detection" , one must also satisfy the criteria that the 
perturbation orbit is well constrained. The exclusion criteria are more robust to such possibilities. 

4.3. Extension to Inclined Orbits 

In order for the exclusion limits to be more broadly useful, it is important to allow for the possibility of inclined orbits 
for the tertiary companions. The present numerical analysis is limited to face-on orbits only due to computational 
limitations; including a search over a range of orbital inclinations introduces two new freely varying parameters, 
increasing the computational time required by roughly two orders of magnitude. Alternatively, a procedure can be 
developed to relate the face-on limits to those for inclined systems. 

For a given face-on system of semimajor axis ap, one can determine the size ai of the inclined orbit model which 
results in the same \ 2 metric; see figure [3] Two cases are evaluated: first, when the observational uncertainty ellipses 
are much larger in one dimension than the other (case for model 71), and second, when the uncertainty ellipses are 
circular (case for model 12). These extreme cases bound the values appropriate for PHASES data. In both cases, it 
is assumed that no perturbation exists and the "wide" binary orbit has been removed; the residuals are centered at 
(x, y) — (0, 0). The \ 2 sum is divided into two terms; one summing along the axis parallel to the larger dimension of 
the apparent orbit of the inclined perturbation model (x), and the other along the smaller dimension (y): 

2_^ {xmi))-xi? | y {y^m-y l ) 2 (3) 
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TABLE 1 

Values of 8q, n 



Q N S Q:N 



Model A 
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25.26 
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Fig. 3. — Inclined tertiary component (perturbation) orbits (II and 12) that fit observations equally well as a face-on orbit model (F). 
Model II corresponds to the case of measurements with uncertainties in dimension x much larger than y\ model 12 corresponds to those 
with equal circular uncertainty ellipses. 



where X((f>(t)) = a cost/) and Y(<p(t)) — acosisin</> (i is the inclination of the orbit) are the perturbation orbit model 
values at time t in each dimension (</> = 2irt / P + tfi ; <f> is a phase offset), and a 2 , and a 2 , are the respective uncertainties 
for each dimension. There are noi-j/ cross terms by choice of basis, explained below. Conditions for which the Xf 
of the face-on orbit equals Xn or X72 °f the inclined models are sought. 

In the case of non-circular uncertainty ellipses, one works in a basis for which the uncertainty ellipse major axis is 
oriented along the x axis. It is assumed the uncertainty ellipses share a single such basis (their orientations do not 
vary) such that a \ 2 minimizing algorithm will align the larger dimension of the apparent perturbation orbit with the 
large uncertainty axis (this also assumes there are no external limits placed on the longitude of the ascending node) . 
For PHASES data, the orientations of the uncertainty ellipses vary slightly from night to night (depending on the hour 
angle of observations); this serves only to lessen the resulting conversion factor, and the present derivation is for the 
more extreme case. In the limit of large uncertainties in the x direction, | X(t) |< an <C cr 2 and the first term in eq. [31 
is unchanged and cancels upon setting Xf ~ Xjx-> leaving 



\p (Yf (U) - yi) = (£n (tj) ~ Vi) 
i=i y i=i v 



(4) 



By geometry, Yp((f>) = ap-sin^ and Yn(<j)) = aji cos i sin <fi. By setting ap = an cosi, \f = Xiv Thus, inclined orbits 
with semimajor axis 1 / cos i times larger will fit the observations equally well as a face-on model. 



(> 



For circular uncertainty ellipses, c 2 = cr 2 



(X f (0(t,))-:r») 2 

t Zji=l c-2 — Z^j=l o-2 

, v-£>/ 2 (y J2 (0(t,))-a,) 2 

X^l ( a -F cos 4>{U) - Xi) 2 
+ Ejil ( a -F sin <t> (**) _ J/*) 2 = Sfc/l 2 ( a -f2 COS (tj) - a;,) 2 

- / 2 f„„ air, ^ _ „,^ 2 



(5) 



+ Z)i=i (a/2cosisin0(ii) - yi)^ 

Assuming the measurements are uniformly distributed in 0, and Xj and j/i have statistically similar distributions p 
which are even functions (such as Gaussians), this relationship becomes 

I^°oo lo* ( flF cos 4> ~ x) 2 p{x)d(j>dx 

+ J™ 00 Jo*( aFsin( t'~yfp(y) d <l ,d y = /^/^(a^cos^-xfp^d^dx 

+ f-oo IT ( a/2 sin <t> cosi - v) 2 P(y)d<f>dy 

f^°ao fo* cos2 ~ 2aF cos cf>x + x 2 ) p(x)d(f>dx 

I-oo la" ( a F sin2 ^ _ 2ap sin ^ + P(y) d ^ = I-oo Jo* ( a ?2 cos 2 - 2a/2 cos 0a; + x 2 ) p{x)d<l>dx (6) 

+ / Q 2lr (a 2 2 sin 2 ^ cos 2 i — 2a/2 sin cos iy + j/ 2 ) p(y)d(/>dy 

J-oo la* a F cos2 "M^d^da; 

+ I-oo la* a F sin2 <h>{y)dct>dy = f^f 2 " a 2 I2 cos 2 (f>p(x)d(f>dx 

+ I-oo It* fl2 2 sin 2 ^ cos 2 ip(y)dcf> 

where in the last step, the second terms (linear in x and y) integrate to zero because p is even, and the third terms 
cancel on the right and left sides of the equation. Finally, the remaining x and y integrals are common to all terms 
and cancel, and the <fi integrals evaluate to identical values for all terms, leaving 

2a% = a% (1 + cos 2 i) . (7) 



Thus, by setting a F = a I2 \/(l + cos 2 i) /2, Xf = X/2- 

Because orbital size is proportional to companion mass, these multiplicative factors can be directly applied to the 
exclusion limits for astrometric companions. To convert between the results for face-on and those including inclined 
models, one need multiply only by a factor between ~ ^/2/(l + cos 2 i) and ~ 1/cosi. The exact value of this factor 
depends on the distributions of observations in time and orientations and aspect ratios of the uncertainty ellipses. 

5. 5 EQUULEI 

PHASES observations of the nearby binary system 5 Equulei (7 Egu.HR 8123, HD 202275, ADS 14773), including an 
updated visual orbit, have been presented bv lMutersnaugh et alJ l|2005j) . The 27 observations reported there have been 
combined with 11 new measurements from the 2006 observing season to search for evidence of astrometric companions 
around either star. Figure0|shows the regions in companion mass-period phase space for which the measurements are 
inconsistent with a perturbation caused by a hypothetical companion in a face-on circular orbit. Companions with 
masses greater than the lines plotted are inconsistent with the PHASES observations; the exclusion regions are those 
above the lines. 

5.1. Application of Detection Limits Procedure 

Implementation of the criteria described by equations ^ an d HI can be demonstrated on this system. The binary 
orbit is short enough that it is well constrained by the PHASES observations alone (Pi < 6 years); the elements of 
the visual orbit are not strongly covariant. Thus, all seven parameters of a visual orbit (P\, e%, i\, u>i, T a %, f2i, a\) are 
allowed to be varied in Keplerian model (Qa = 7). The resulting Xa is 1090.83. The 38 PHASES measurements are 
each two-dimensional, the expected x\ is thus Za = 2 x 38 — 7 = 69; under the hypothesis that the single-Keplerian 
model is correct, \\ is m excess by a multiplicative factor of 1090.83/69 = 15.81. The measurement uncertainties 
themselves are too small by a factor of -\/ 15.81 = 3.98, and the median minor-axis formal uncertainty of 7.7/Ltas is 
corrected to 30.6/ias. 

Model B is a double-Keplerian orbit given by superposition as 

Vohl = tK-b (Pi,ei,ii,wi,T 0) i,f2i,ai) + 

rBa-Bb.C.O.L! (-P2, e2, fa, 0J2, T ,2, ^2, a 2 ) ■ (8) 

Here, the A-B orbit is the slow, "wide" orbit and the Ba-Bb orbit is the "narrow" perturbation orbit center-of-light 
motion. When converting a 2 to M-gb, assumed values for Me a = 1.19Mq and the system distance d — 18.386 pc are 
used, and it is assumed that component Bb is faint (i.e. the Ba-Bb center-of-light is located at Ba, and a 2 = aB a , 
corresponding to just the reflex motion of the star). In S Equ, both stellar components A and B mass roughly 1.19Mq; 
the model is symmetric for companions to A or B and the derived limits on companion masses are identical for the 
two — one search eliminates companions around both stars, in this case with nearly equal mass limits. 
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To evaluate whether model B (equation[5J) better fits the observations, the x 2 metric is used to evaluate the goodness- 
of-fit at each point in a grid of hypothetical companion masses and orbital periods. The grid was stepped logarithmically 
in companion mass from 10 _1 to 10 3 Mj with step size log (MBb/-Wj) = 0.1. The companion period was sampled 
at values of Pi = 2/T/K, where / is an excess factor (here 3) for finer sampling, T is the span of the observations 
(~ 760 days), and K is a natural number up to that for which Pi = 3 days. Pi was additionally evaluated at values 
of K= 1/2, 1/3, ...,1/9. 

At each grid point, the seven parameters of the wide orbit are seeded with those from the best-fit single-Keplerian 
(i.e. model A); the perturbation orbit's parameters are seeded in period and (center-of-light) semimajor axis as deter- 
mined by the phase-space grid. Additionally, the perturbation orbit is seeded in several values of tti (separated by 40 
degrees starting at and increasing up to 320), and at inclinations of either zero or 180 degrees to allow for face-on 
orbits with clockwise or counter-clockwise motion; all combinations of these £li and i% seedings are explored. The 
values of the other secondary orbit parameters are seeded at zero. 

The model is allowed to relax from each set of initial seedings using a downhill \ 2 fitting algorithm in which 
the seven parameters of the "wide" orbit and in the "narrow" orbit are allowed to vary. Note that the narrow 
system inclination is not allowed to vary — this restricts the search to including only face-on orbits, and the inclination 
parameter does not add to the count of "free parameters" in the model. While the period and semimajor axis (and 
thus mass) of the perturbation orbit are not allowed to vary in the downhill fitting process, these are varied across the 
grid and are to be counted in the total number of model parameters. Note that with eccentricity held at zero (circular 
orbits) and inclinations restricted to zero and 180 degrees, fli is degenerate with T Q _i and u>i, thus these need not be 
varied. The total number of degrees of freedom in the double-Keplerian model is thus Qg = 10; seven are from the 
wide orbit, plus the period P%, orbit phase f^, and semimajor axis ai (or companion mass) of the perturbing orbit. 
From this, = Zj, — 3 = 66. 

The Na exclusion regions are defined as those for which the double-Keplerian model x% ^ s greater than 1090.83 x 
(66+<5io,jv) /69. The Na detection regions are calculated from equation|21as being those for which the double-Keplerian 
X% is less than 1090.83 x 66/(69 + 67, n). The la (x% = 1226) and 3a (xb = 1469) exclusion contours in mass-period 
phase space are plotted in figure^ also shown is the la detection contours at Xb = 933; the smallest value of Xb ls 
883, well above the 2a detection criteria of Xb — 864. It is concluded these detection regions are results of statistical 
happenstance. 

5.2. Comparison to Previous Work 

The same procedure has been applied to the non-PHASES astrometric measurements reviewed bv lMuterspaueh et alJ 
1)2005(1 : see this work for a discussion of the relative weightings of these measurements. These 170 measurements cover 
a much longer time span (of a century) than the PHASES observations, but at lower precisions. In this case, the value 
of T used to determine the companion periods sampled was chosen to be half the period of the wide orbit found in 
the best fit single-Keplerian model (Pi/2 = T ~ 1042) days. The best fit seven-parameter single-Keplerian model 
results in x\ = 426.5 and = 340 — 7 = 333 degrees of freedom. Thus, the la and 3a exclusion level for Xb of the 
ten-parameter double-Keplerian model are 437.5 and 457.2, respectively. For Pi < Pi, no values of x\ m the grid are 
less than 400, which is between the 2<r and 3a detection thresholds (this lowest point is near Pi ~ 10.25 days and 
M B b ~ 800Mj, which is strongly disfavored by the PHASES data). 

5.3. Inclined Orbits 

At this time, it is computationally prohibitive to search in additional parameters (for example, to add an array of 
seeds in inclination and allow this parameter to also vary, or to search for eccentric orbits). It is noted that there is 
always an orbit for which astrometry is completely insensitive to companions of any mass or period — high eccentricity, 
edge-on orbits pointing directly at the observer (i.e. periastron passage within the line of sight); combination fitting 
with radial velocity observations lifts this problem and can be useful when contemplating universal exclusions of 
perturbing companions. Evaluating the phase-space grid described above using just the PHASES astrometry and 
a 10 parameter model requires roughly 3 weeks on a modern processor; each additional parameter is expected to 
increase this by roughly an order of magnitude. A need is identified for an alternative approach than the "brute-force" 
algorithm described here, in combination with an evaluation of its performance and reliability in comparison to the 
method described here. 

To test the analytical conversion factor between detection limits for face-on and inclined systems of M4.3I the PHASES 
observations of S Equulei were reanalyzed allowing inclinations to vary within 45 degrees of face-on (both prograde 
and retrograde orbits were allowed) over a very limited subset of perturbation periodicities. The predicted conversion 
factor is between i/4/3 and a/2; in practice, a factor of ~ 1.35 — 1.4 is appropriate for the PHASES observations. 

5.4. Fit Residuals 

Figure [5] shows the continuous (integrated) distribution function of the residuals to the 2-body Keplerian fit for the 5 
Equulei system. For each measurement, the residuals were measured along directions parallel to the formal uncertainty 
error ellipse major and minor axis in order that the uncertainties are not covariant. Each were then normalized by 
that measurement's formal uncertainty estimate. The distribution of this normalized set of residuals can now be 
considered to determine whether a multiplicative scale factor is an appropriate assumption for handling the observed 
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Fig. 4. — The S Equulei Mass-Period companion phase space. Companions in the regions above the plotted exclusion curves with face-on 
orbits are not consistent with the observations. Companions as small as two Jupiter masses can be ruled out by PHASES observations. 
The most significant detections are less than 2<r significant, and are not likely astrophysical in origin. 



excess scatter, in which case the residuals should have a Gaussian distribution. The continuous distribution function 
is fit to the integral of a Gaussian: 

1 + erf (x/V2a) 

2 (9) 

where a is a measure of the excess scatter over the formal uncerta inties (the best fi t value for a is 2.94 ± 0.01). A 
one-sided Kolmogorov-Smirnov test (see for example section 14.3 of iPress et alJ il992') shows the residual distribution 
agrees with Gaussian statistics with 95% likelihood. The distributions match well and the assumption of Gaussian 
errors and multiplicative uncertainty scale factors is valid. Treated as separate sets, the scale factors for the major-axis 
and minor-axis residuals are similar. 




Normalized Residual 

Fig. 5. — The distribution of normalized residuals for a 2-body Keplerian fit to S Equulei. Also shown is the theoretical distribution for 
Gaussian noise. 

6. KNOWN TRIPLE STARS 

This procedure is tested on two binaries that each are known to host tertiary c ompanions, k Pegasi and V81 9 Herculis. 
The PHASES astrometric orbits of these triples have been previously examined ( Muterspaueh et alJ2006blal) : the phase 
space exploration algorithm detects these perturbations and their statistical significances are now presented. 



6.1. k Pegasi 

PHASES observations of the well-known triple system k Pegasi (10 Peg, HR 8315, HD 206901, ADS 15281), 
from which the orbits of both the long and short period subsystems were determined, have been presented by 
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iMuterspaueh et aP (|2006b). The 52 PHASES measurements presented there are combined with three measurements 
from 2005 creating a set spanning ~ 1046 days. In the previous investigation, both a noise floor and multiplicative 
factor were explored for explaining excess systematic noise; the present investigation evaluates only the latter case. 

A mass-period grid was constructed similar to that in the case of 5 Equulei, differing in the following ways: the P2 
oversample factor / is 1 rather than 3, the companion masses explored ranged from 10° to 10 4 Mj, and the step size 
in companion mass was log (M-Qh/Mj) = 0.2 rather than 0.1. The 424 non-PHASES astrometry data tabulated in the 
previous investigation were similarly evaluated over the same grid; because these span a time much longer than the 
wide orbit period, T was set to half the best-fit wide system period P\ (for evaluation of the PHASES data, T is the 
timespan they cover). For both the PHASES and non-PHASES data sets, the wide orbit is well-constrained without 
strong covariances. Thus, for both cases Qa = 7 and Qb — 10. The companion Bb mass to 02 conversion assumed 
a stellar mass Mb s — 1.662M Q and distance to the system of 34.60 pc, as determined by the previous investigation. 
Note that this analysis is only to confirm the detection of the known 5.97-day period companion Bb; at this time no 
attempt is made to search for fourth components to the system, which may require much more complicated modeling 
than the simple superpositioning of independent Keplerians that has been used here. 

In the PHASES evaluation, the smallest value of Xb ^ s f° un d at P2 = 5.978 days; in this region the step size is 
AP 2 ~ 0.017. This is consistent with the best-fit value of P 2 of 5.9714971 ± 1.3 x 10~ 6 days. The depth of this 
detection corresponds to a detection at the 31.4<r level. Conversely, the smallest value of Xg for the non-PHASES data 
is at P2 — 19.14 days with depth not quite reaching the 3a detection level. The smallest value with 5.9 < P2 < 6.1 
days is just barely smaller than the la threshold (near 1.1a). Figure El plots the relevant phase-space contours. 

Upon removing the face-on, circular orbit criteria for the narrow orbit (seeding a \ 2 downhill fit with the best set 
of parameters found in the phase-space analysis), the fit is considerably improved, Qg — 14, and the second orbit 
improves the fit at the 106cr level. 

6.2. V819 Herculis 

PHASES observations have previously deter mined the astrometric orbi t of the short-period subsystem in the V819 
Herculis (HR 6469, HD 157482) triple system ijMuterspaueh et al.ll2006a(l . The narrow (Ba-Bb) pair of the V819 Her 
system is oriented edge-on and show eclipses. In that analysis, 31 measurements were reported (of which six were not 
used in model fitting as they were taken during eclipse of the narrow pair) ; one additional measurement from summer 
2005 is added to this set. In the present analysis, measurements taken during eclipse are not omitted to simulate a 
blind search. 

When only the PHASES data are fit, the span of the measurements (T ~ 476 days) is much shorter than the period 
of the wide orbit (Pi ~ 2020 days). Thus, Pi, e\, and a\ are strongly covariant; by fixing just one of these (here, 
Pi = 2019.787 days) at a nominal value supported by non-PHASES observations, these covariances are lifted and the 
fitting algorithm converges much more rapidly. Because that parameter is not optimized for each fit, it does not count 
to the number of free parameters; for the PHASES analysis, Qa = 6 and Qg — 9; this is not necessary when the 
non-PHASES observations are included, for which Q_a = 7 and Qb = 10 as normal. 

As in the case of k Pegasi, the mass-period phase space grid was modified by setting the P2 oversample factor / 
to 1 rather than 3, the companion masses explored ranged from 10° to 10 Mj, and the step size in companion mass 
was log (M-Bh/Mj) = 0.2 rather than 0.1. Additionally, the smallest value of P2 sought was 1 day rather than 3. The 
distance is assumed to be 67.96 pc and Mb s = 1.43M©. Again, the aim of the search is only to confirm detection of 
the known companion Bb, and no attempt is made to find fourth components at this time. 

Because the algorithm has been developed to specialize in face-on orbits, it is not optimally suited to analysis of 
the V819 Her system, for which the short-period perturbation system is observed to eclipse. However, the strength 
of the algorithm as a "wobble-detector" is demonstrated by applying it to this system for which it is quite non- 
optimally designed. The two deepest dips in the x 2 function occur at periodicities of 2.2284 days (AP2 ~ 0.0052 
days; 10.5a) and 1.8056 days (AP2 ~ 0.0034 days; 10.1a). These are the only two points at the 10a level; the 1.8056 
day periodicity appears to be only an aliasing; the RV and eclipse photometry coincide with the (more significant) 
2.2284 day periodicity (the best fit period to all data is P2 = 2.2296337 ± 1.9 x 10~ 6 days, in agreement with the 
2.2284 periodicity, to within the analysis step size). This short-period, 110/uas perturbation is readily detected, but 
is far beyond the ability of other astrometric measurements; with this periodicity identified, one can better refine the 
fit while allowing other parameters to vary, obtaining the correct orbital configuration for the short-period pair and 
improving the detection to the 14.7a level. 

The 34 speckle interferometry data tabulated in the previous paper are evaluated separately (in this case Pi is 
allowed to vary freely) . As shown in figure the speckle observations are approximately two orders of magnitude 
away from detecting the short-period perturbation. Interestingly, the llO^as amplitude perturbation would correspond 
to a planetary mass companion if only its period were several months rather than several days. 

7. CORONAE BOREALIS 

(3 Coronae Borealis ("Peculiar Rosette Stone", 3 CrB, HR 5747, HD 137909) is a binary whose primary is of spectral 
class F0p, part of the family of peculiar A stars with strong magnetic fields. (3 CrB is often grouped with 7 Equulei 
and a 2 CVn as prototypes of this class. Its magnetic field has been extensively via measurements of Zeeman line 
s plitting, and the inclination of its rotation axis and the angular offset of its magnetic field have been measured. 

iNeubauerl l)1944f) studied the binary using RV measure ments and found evid ence for a third body with orbiting the 
primary with a period of nearly a year (P2 ~ 320 days). iKamper et alJ l)1990j) obtained updated RV data which did 
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Fig. 6. — k Pegasi Mass-Period Companion Phase Space. 
(Top) The ltr exclusion contour based on previous differential astrometry data, with the 20cr PHASES detection contour. Note that only 
at the 5.97 day companion period does the 20a PHASES contour exist. 

(Bottom) The 10-, 20-, and 30-<r detection contours for PHASES observations of the n Pegasi triple star system are plotted for the limiting 
case of face-on Ba-Bb subsystem orbits. While the Ba-Bb subsystem is not in fact face-on, the perturbation is detected with the same 
periodicity as that found by RV observations. Once this periodicity has been identified, the face-on, circular orbit constraint can be 
removed, and the detection becomes 106cr significant. Also plotted is the lcr exclusion contour based on previous differential astrometry 
data. A "detection" at the ltr level in the previous differential astrometry is seen in the plot; however, this is likely due only to statistical 
happenstance. 



not show this perturbation. However, their analysis showed the periodicity in lNeubauerf s data is statistically different 
than one year. They suggested the perturbation was real rather than instrumental, and proposed that th is orbit was 
previo usly nearly face-on, somehow re-orienting itself to be perfectly face-on between the two epochs. ISoderhielml 
(1999) concluded from Hipparcos astrometric observations that such a proposed companion cannot exist. 

Forty-two PHASES measurements have been collected over a span of ~ 804 days. Additionally, 102 non-PHASES 
differential astrometry measurements have been identified (when these are analyzed, T is set to half the best fit value 
of Pi in the 2-body single-Keplerian case). As was the case for V819 Her, Pi, e\, and a\ are strongly covariant when 
only PHASES data are considered (T ~ 804 days versus Pi ~ 10.27 yr). Again , Pi is fixed at a value that is not based 
on the PHASES observations; Pi = 10.27 yr is adopted from iTokovininl l|1984D . and Qa — 6, Qb = 9 (7 and 10 when 
non-PHASES data are included). 

The mass-period grid is set up in the same manner as for 8 Equulei. Based on analysis of Hipparcos observations 
bv ISoderhielml l|1999f) . the distance is assumed to be 34.12 pc (n — 29.31 ± 0.82 mas) and the component masses are 
1.77 ± 0.16 and 1.21 ± 0.11M Q . The conversion of tertiary mass to semi-major axis assumes a stellar mass of 1.21M Q . 
Radial velocity measurements are available only for one component of the binary; this prevents one from determining 
the masses and distance from the orbit alone. 

The excess multiplicative factor is Xa/^a = 9.9 ~ 3.2 2 for the PHASES observations. The lcr and 3c exclusion 
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1o Exclusion, Previous Observations 

10a Detection, PHASES 

Fig. 7. — The V819 Her Mass-Period companion phase space shows a 10a significant detection in the PHASES data that is far beyond 
the detection threshold for previous measurements. Although the tertiary companion's orbit is in fact edge-on, the face-on code detects 
the wobble at the correct period. After removing the face-on orbit restraint, the double-Keplerian model is a 14.7cr improvement over the 
single-Keplerian model. 



contours for both the PHASES and non-PHASES data sets are shown in figure |S| In both cases, the Pi ~ 320 day 
perturbation is excluded; for the PHASES observations, a companion of mass as small as 4Mj are excluded. The 
mass limits should be increased by ~ 1.77/1.21 ~ 1.46 when considering the primary star as the companion host (this 
estimate assumes the companion is an insignificant p art of the total mass, reasonable in the case of the lower limit 
exclusion region). It appears that the perturbation in iNeubauerf s data set is not astrophysical in origin. 

There are regions of the mass-period phase space for which perturbations from a tertiary companion improve the 
fit by 3<T or more. The smallest value of x% corresponds to a 5.2er detection at Pi ~ 20.9 days. When the restriction 
to face-on, the narrow orbit's parameters are strongly covariant and poorly constrained: for example, the formal 
uncertainty on the eccentricity is 0.5, implying it is not constrained at all, and the Eulerian angles all have formal 
uncertainties greater than 360 degrees. While the double-Keplerian model can formally improve the fit compared to 
the single version, the double-Keplerian model is too detailed to be properly constrained by the current observations. 
In order to claim a true detection, multiple criteria must be met: the double-Keplerian model must better fit the 
data, and the perturbation orbit must be well-constrained. Thus, it is possible that the "detection" periodicities are 
astrophysical in nature, but the data sampling in the present set prevents conclusive study, and aliasing causes multiple 
periodicities to exist. 
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Fig. 8. — The mass-period phase space for tertiary companions to /3 Coronae Borealis shows that a proposed massive object with 
period ~ 320 days does not exist. Some mass-period combinations do result in double-Keplerian fits that are an improvement over the 
single-Keplerian model, but the inner-body's orbit is poorly constrained and it is not clear whether these are astrophysical in origin. 



8. 13 PEGASI 



13 Pegasi (HR 8344, HD 207652) has been observed 25 times over a span of T ~ 441 day s in the PHASES program. 
The span is short compared to the binary period (Pi ~ 26.132 yr; lHartkopf et al.lll989|) and again Pi, e%, and a\ 
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are covariant; in the analysis Pi is thus held fixed at this nominal value (Qa = 6, Qb = 9). So derhielml l)1999[) used 
Hipparcos measurements to determine the system parallax of 29.6 ± 0.8 mas, a total mass of 2.67Mq, and photometric 
mass ratio Mb/Ma = 0.80. Computations assume a distance of 33.8 pc and the tertiary companion orbits the secondary 
at mass 1.19M Q (the results can be scaled appropriately for the primary at mass 1.48M . 

The best fit Xa = 148.8 with Za = 44 degrees of freedom; assuming model A is correct, the multiplicative excess 
factor is 3.38 ~ 1.84 2 . The median formal minor axis uncertainty of 10.2/zas is thus adjusted to 18.8/xas repeatability 
from night to night. The mass-period grid is set up as in the case of 8 Equulei. Figure [5] shows that massive planets 
in 20 day to 3 year period face-on circular orbits would perturb this binary by more than the observed scatter in the 
PHASES data. Planets as small as two Jupiter masses are ruled out in ~ 4 month period orbits. A number of regions 
show lcr consistent detections, the best of these is at only the 1.7a detection level. 
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Fig. 9. — The 13 Pegasi Mass-Period companion phase space shows PHASES observations can rule out tertiary objects as small as two 
Jupiter masses. A few mass-period combinations introduce slight improvements over the single-Keplerian model, but none of these are 
more significant than 1.7<x, and are probably not astrophysical in origin. 



9. 20 PERSEI 

The binary 20 Persei (HR 855, HD 17904, ADS 2200) has been observed 26 times by PHASES. At K = 4.3 it is 
one of the fainter PHASES targets. The observations span T ~ 877 days, a small amount compared to the orbital 
period of Pi ~ 31.528 yr. Covariances between Pi, e\, and a\ are extremely strong. Using values from lDocobo et all 
(2001), both Pi (31.528 yr) and ei (0.753) are fixed {Qa — 5, Qb — 8). The system distance is assumed to be 72.1 pc 
from the Hipparcos parallax; combining this with the orbital period and semimajor axis an average component mass 
of 2.1Mq is computed and these values are used in converting tertiary mass to 02. The mass-period grid is set up 
following the same procedure as for 8 Equulei. For 20 Per, Xa/^a = 6.57 ~ 2.56 2 is a relatively small value. However, 
the average formal minor-axis uncertainty is 21/ias; coupling this with the large distance and masses implies that the 
tertiary detection limits are not particularly sensitive. 

The PHASES measurements rule out brown dwarf mass companions in a variety of orbital periods, as shown in 
figure El Detections are found for several different companion periods at the la detection level. The smallest value 
of x% appears at the 2.9a level; as in the case of j3 CrB, when t he face-on and c ircula r restrictions are lifted from 
the perturbation orbit, its parameters become poorly constrained. lAbt and Levvl {1976) proposed the existence of a 
tertiary companion with period Pi — 1269 ± 70 days; this periodicity is very close to the steep cutoff in the phase space 
contour due to the finite span of PHASES observations; at this time no conclusions are made about the existence of 
such a companion. 

10. rj CORONAE BOREALIS 

Seventeen PHASES observations of 77 Coronae Borealis (2 CrB, HR 5727, HD 137107, ADS 9617) have been collected 
over a span of T ~ 680 d ays. As with 20 Pe r, Pi and ei are fixed at values of ~ 41.585 yr and 0.2620 determined 
by spe ckle interferometry ((Mason et al.lll999T) to lift strong covariances with a%. From Hipparcos data. IS5derhielm1 
( 1999) determined a parallax of 53.5 ± 0.9 mas, total mass 2.41M Q , and photometric mass ratio of 0.96. From these, 
the tertiary orbit to 02 conversion is made assuming a distance of 18.69 pc and stellar mass of I.I8M0 (the primary's 
mass is roughly 1.23Mq and the tertiary companion limits are similar for primary and secondary). The mass-period 
grid is set up as in the case of 8 Equulei. 

Companions as small as 3-4Afj in orbits of ^100-700 days are inconsistent with the PHASES observations. As in 
the case of (3 CrB, detection contours deeper than 3a appear (here at periods less than 10 days). The deepest of these 
is at the 4.8a level. Similarly to (3 CrB, the data sampling does not allow all parameters of the perturbation orbit to 
be well constrained when the restrictions of face-on and circular orbits are lifted (in particular, ei is not constrained 
at all). This may be astrophysical in origin, but the current data cannot properly address the underlying orbit. 
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Fig. 10. — A range of brown dwarfs companions can be ruled out by the 20 Persei mass-period phas e space. Unfortunatel y, the PHASES 
observations do not yet have long enough time coverage to confirm the 1269 day period companion of Abt and Levy ( 197G). Owing to the 
far distance of the system and higher mass stars, the PHASES observations are not particularly sensitive to planetary companions in this 
system. 
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Fig. 11. — The r\ Coronac Borealis Mass-Period companion phase space shows planets as small as three Jupiter masses are inconsistent 
with PHASES observations. A nearly 5<r significant detection is found at short periods, but observational coverage does not yet allow the 
orbit to be strongly constrained; it is unclear whether this is merely a statistical fluke, perhaps related to the data sampling function (both 
dips are near 7 days, a natural sampling time for PHASES observations). 



11. HR 7162 

HR 7162 (HD 176051, ADS 11871) is a GOV primary with K1V secondary. This color difference appears to affect the 
PHASES data through the differential dispersion mechanism because the excess scatter Xa/%A 1S much larger than 
for other stars (~ 88 ~ 9.4 2 ) . As with 20 Per and r\ CrB, both P\ and ei are fixed at values of ~ 61.15 yr and 0.258 
determined bv iHeintl l)1994f) to lift strong covariances with a%. Over a span of ~ 773 days, 32 PHASES observations 
have been collected. The tertiary mass to conversion ass umes a distance of 15 parsecs and star mass 0.71M Q . These 
quantities are derived from the Hipparcos-based solution of lSoderhielml (l999), who found a parallax of 66.7 ±0.6 mas, 
mass sum 1.78-Mq and photometric mass ratio 0.67 (the primary mass corresponds to 1.07-Mq). The mass-period grid 
is constructed as in the case of 5 Equulei. The smallest values of x| correspond to a 2.7c detection. Companions as 
small as 2Mj are inconsistent with the PHASES data. 

The physical properties of the HR 7162 binary (low mass stars, close to Earth) make it an exciting candidate for an 
exoplanet search as very low mass companions can be detected for a given astrometric performance levels. Unfortu- 
nately, the PHASES excess scatter factor for this binary is currently much larger than that typical in other systems 
studied, and appears to be random noise. This is likely due to the large color difference of the binary components, 
potentially leading to differential dispersion issues. It is anticipated that the current engineering improvements will 
benefit this system in particular, in which case it will be an extremely exciting study for the exoplanet search, with 
the ability to detect (or rule out) planets as small as a fifth of a Jupiter mass orbiting either star. 

12. CONCLUSIONS 

PHASES observations are able to exclude tertiary companions with masses as small as a few Jupiter masses in 
several binary systems. The phase-space searching algorithm correctly identifies the natures of two previously studied 
triple star systems and finds the correct periods for the short period pairs. While the number of systems that can 
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Fig. 12. — Contours in the companion mass-period phase space for HR 7162 show that despite high noise levels from systematic processes, 
companions as small as two Jupiter masses can be excluded in face-on orbits. After engineering improvements to remove these systematic 
noise sources, constraints will be placed on sub-Jupiter massed planets by second stage PHASES observations. 



currently be examined is too small to make strong conclusions about the frequency of giant planets in close binary 
systems, by the end of the PHASES program enough measurements will be made on all target systems to address this 
important question. Systematic effects currently limit the observed precision; if overcome, companions as small as a 
fifth of a Jupiter mass might be detected in favorable systems such as HR 7162. 
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